function RK=rk4gral(a,b,ya,M,func)

h=(b-a)/M;
T=zeros(1,M+1);
Y=zeros(1,M+1);
T=a:h:b;
Y(1)=ya;

for j=1:M   
	K1=h*feval(func,T(j),Y(j));
	K2=h*feval(func,T(j)+(h/2),Y(j)+(K1/2));
	K3=h*feval(func,T(j)+(h/2),Y(j)+(K2/2));
	K4=h*feval(func,T(j)+h,Y(j)+K3);

	Y(j+1)=Y(j)+(1/6)*(K1+2*K2+2*K3+K4);
end

RK=[T' Y'];
